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Abstract 

We examine the problem of the collapse and fragmentation of molecular clouds with 
a Gaussian density distribution with high resolution, double precision numerical 
simulations using the GADGET-2 code. To describe the thermodynamic properties 
of the cloud during the collapse -to mimic the rise of temperature predicted by radia- 
tive transfer- we use a barotropic equation of state that introduces a critical density 
to separate the isothermal and adiabatic regimes. We discuss the effects of this cri- 
tical density in the formation of multiple systems. We confirm the tendency found 
for Plummer and Gaussian models that if the collapse changes from isothermal to 
adiabatic at earlier times that occurs for the models with a lower critical density, 
the collapse is slowed down, and this enhances the fragments' change to survive. 
However, this effect happens up to a threshold density below which single systems 
tend to form. On the other hand, by setting a bigger initial perturbation amplitude, 
the collapse is faster and in some cases a final single object is formed. 



1 Introduction 

The protostellar objects that begin their main sequence and pre-main sequence 
are mainly distributed in binary and multiple systems that suggests they were 
formed during the process of collapse and fragmentatio n of m olecular clouds with 
dense cores and gas envelopes, see i Sigalotti (2Q0li;[ Tohline (2002.) and refer- 



ences therein. In recent years several authors have considered different realiza- 
tions of molecular clouds to study their collapse and fragmentation. Although most 
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fragmentation calculations apply to initially uniform conditions, see for instance 
fBodenheimer ( 200 0), it is clear from the obse r vations that mo l ecular c l oud co res are 
centrally condensed rWard-Thompsonl(l2002b : l Andrel(ll998l) : l Mottd ( Il998l) . Thus, 
a number of collapse models starting from centrally condensed, Gaussian density 
profiles have also been made. A particular co mputa t ionally demanding isothermal, 
Gaussian cloud model was first calculated in I BossI (119911) . and thereafter recalcu- 
lated by other authors as a further test case to check both the likelihood of fragmen- 
tation during the isothe r mal collapse phas e and the reliability o f the numer i cal cod e 
results I Burkertl (ll996l): I Truelovel (Il997l) : I BossI (Il998l |2000|) : I Sigalottil (l2001bl) : 



SigalottilOOOld) . 



So far the great majority of this research has concentrated upon the early phases 
of star formation, when the collapse is dynamical, first isothermal and then non- 
isothermal. However, precise knowledge of the dependence of temperature on den- 
sity at the transition from isothermal to nonisothermal collapse requires solving the 
radiative transfer problem coupled to a fully self-consistent energy equation. How- 
ever, the full non-isothermal computation represents a severe computational burden 
imposed by solving the radiative transfer equations at high spatial resolution, even 
in the Eddington approximation. T herefore, it ha s been common to use instead a 
simple barotropic equation of state I Boss that clearly simplifies the com- 

putational problem and it turns out to be a good approximati o n for the dynam ical 
collapse of the molecular clouds, see I Arreaga-Garcial (120071) : I Arreagal (120081) . In 
the present work, a barotropic equation of state is assumed to simulate the transition 
from isothermal to adiabatic collapse. The motivation of this study is to investigate 
the sensitivity of fragmentation to the effects of thermal retardation by varying the 
value of th e critical dens ity at whi c h non isothermal heating is assumed to begin. In 



Arreaga-G arcfa (2007 J; I Arreagal (12008b it is studied the evolution of a Gaussian 



density profile, and found that by diminishing the critical densit y it enhances the 
fragm entation. A similar result was found for Plummer models in I Arreaga-Garcia 

(1201 Ol) . 

The p r esent work is a continuation of the analysis in I Arreaga-Garcfal (l2007b : 



Arreagal (|2008|) in which we employ a Gaussian density profile and perform 



the same type of numerical computations but now using double precision in the 
GADGET-2 code. Particularly, we study the effect of varying the critical density of 
the barotropic equation of state in the collapse and fragmentation of the molecular 
protostar. We also analyze the effect to increasing the initial perturbation ampUtude. 



2 Initial conditions and collapse models 

According to astronomical observations, the regions from which stars are formed 
consist basically of molecular hydrogen clouds at a temperature of ^ IQK. There- 
fore, the ideal equation of state is a good approximation to account for the thermo- 
dynamics of the gas in these clouds. The cloud m odels are based on the standard 



isothermal test case, as in the variant considered in l BurkertI (119931) . However, once 
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gravity has produced a substantial contraction of the cloud, the opacity increases, 
the collapse changes from isothermal to adiabatic and the gas begins to heat. To in- 
clude this rise in te mperature int o our calculations, we use the barotropic equation 
of state proposed in I Boss (l2000l) . 

In order to correctly describe the non-isothermal regime, one needs to solve the 
radiative transfer problem coupled to the hydrodynamic equations, including a fully 
self-consistent energy equation to obtain a precise knowledge of the dependence of 
temperature on density. The implementation of radiative transfer has already been 
included in some mesh-based codes. In SPH, the incorporation of radiative transfer 
has in general not been very satisfactory, perhaps with the exception of reference 
fwhitehous e (200©, in which they used the f-limited diffusion approximation to 
model the collapse of molecular cloud cores. These authors suggested that there 
are important differences in the temperature evolution of the cloud when radiative 
transfer is properly taken into account. 

However, after comparing the results of the simulations perf ormed by 

Arreaga-Garcia ( 2007 ): Arreaga ( 2008 1 with those of reference Whitehouse ( 20061) 



for the uniform density cloud, it is concluded that the barotropic equation of state in 
general behaves quite well and that we can capture the essential dynamical behavior 
of the collapse. The simulations in this work are consequently carried out using the 
following barotropic equation of state; 



(1) 



where 7 is the adiabatic exponent in the opacity thick regime and A" is a constant set 
by A" = cj^uPf.^/, where Pcit defines the critical density above which the collapse 
changes from isothermal to adiabatic, and for a molecular hydrogen gas the ratio of 
specific heats is 7 = 5/3, because we only consider translational degrees of freedom. 
With the above prescriptions, the local sound speed becomes 



P 

Pcrit 



y-1 



-,1/2 



(2) 



so that c ~ Ciso when p ^ PcrU and c w Cad = /'''^Cuo- 

The molecular cloud collapse simulations in this work begin with i nitial con- 
ditions in accordance to the thermodynamic model proposed in iBurkert (1993j). 
Accordingly, the models start with a spherical cloud of mass M = \Mq, radius 
R = 4.99 X lO'^cm ~ 0.016pc, and at a temperature T = lOK. The initial model 
is composed by an ideal gas with an average molecular weight ^ 3. We have 
chosen the initial sound speed c,,,, and angular velocity Oo in such a way that for 
all models the initial ratio of thermal and rotational energies to gravitational en- 
ergy are such that a = -Ethemi/I^gravl ^ 0.26 and j3 = -Erot/I^gravl ^ 0.16. The gas 
isothermal sound speed is ctso ^ 1.90 x lO^cms^', and the average free fall time is 
5.10 X 10"s. Additionally, we impose a small perturbation to the density profile of 
the following form; 

p = po[l+acos{m^)], (3) 
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where m is an integer number, the azimuthal angle around the z— axes, and a is 
the perturbation ampUtude. 

The ch osen density profile is Gaussian wi t h the a bove-mentioned initial condi- 
tions, as in . 



Arreaga-Garcial (120071 ): I Sigalottil(l2001al) : 

p(r) =p,.exp 



(4) 



where = 1.7 x 10— 17grcm ^ is the initial central density and b ~ 0.578/? is a 
length chosen such that the density is 20 times smaller there. On the other hand. 



-12,-1 



solid-body rotation is assumed at the rate of COq = 1 .0 x 10 s 



3 Numerical methods 



The computations of this work were performed using the parallel code GADGET- 
2, which is described in full in ISpringel I (120051) . The code is suitable for studying 
isolated, self-gravitating systems with high spatial resolution. The code is based on 
the tree-PM methods for computing the gravitational forces and on standard SPH 
methods for solving the 3D Euler hydrodynamics e quations. For a re view on the 
theory and applications of SPH we refer the reader to I MonaganI (120051) . 

In order to set up the initial particle distribution, we first define a Cartesian box 
with sides equal to twice a specified radius Ri, > R ^ 4.99 x 10^^ cm, and with 
its geometrical center coinciding with the origin {x — y = z = 0) of a Cartesian 
coordinate system. The box is then subdivided into regular cubics cells of volume 

= AxAyAz each. The spherical cloud is then copied within the box by placing an 
SPH particle in the center of each cell at distances d <R from the origin, so that the 
region outside the sphere is empty. A little amount of disorder is added to the regular 
distribution of particles by shifting each particles a distance A/ 4- from its cell-center 
location and along a specified direction, which is chosen randomly among the three 
Cartesian axes. We defined the mass of particle / at location (x, , y,-, z,) to be m, — 
p{xi,yt, zi)A^ , where 



p{xi,yi,Zi) = Pc&W 



xj+yj + zf 



(5) 



Solid-body rotation about the z — axis is assumed in a counter clockwise sense by 
assigning to particle ; an initial velocity given by v; = {(O^Xi, — a)oy,-,0). Finaly, the 
bar mode density perturbation given by equation (|3) applied by modifying the mass 
of particle / according to ot,- ot, [1 + acos (ot0,)], where 0,- denotes the azimuthal 
position of that particle. The computations were performed in the parallel cluster of 
the National Institute of Nuclear Research-Mexico, equipped with 28 AMD Quad- 
core (64 bits) Opteron Barcelona processors. 
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In this section we present the resuhs obtained using the GADGET-2 code of the co- 
llapse and fragmentation of Gaussi a n mol ecular clouds. The coll apse of the Gaus - 
sian cloud first calculated in Bos s later by other authors I Burkerd (1 1996b: 

Truelovel(ll997l) : l Bossl(ll998b :l Sigalotti(2001h'): ISigalottil ( l2001ch : l Arreaga-Garcia 
(l2007h using high spatial resolution , and in the present work we perform the same 



computations as in I Arreaga-Garcial (12007b . but adding double precision and using 



10^ SPH particles in each simulation. 

We present four different cases for the barotropic collapse with the parameters 
as shown in table [T] The results are illustrated with iso-density contour plots for 
a slide at the equatorial plane of the cloud in figures [T]|4] The bar located at the 
bottom of the plots shows the logjo density range at a time t and normalized with 
the initial central density pc- A color scale is then associated with the value of logjg. 
For instance, the color scale uses yellow to indicate higher densities, blue for lower 
densities, and green and orange for intermediate densities. 

The free fall time f/-/ ~ ^37r/ (32Gpf) sets a characteristic time scale for the 
collapse of protostellar clouds which is gi ven in terms of the cent ral density and it 
is the same for all models, see for instance ! Arreaga-Garcia I (l2010l) . 



Table 1 Gaussian collapse models. The model types are explained as follows: The letter G refers to 
Gaussian, the number 6 refers to ten millions particles, the letters A, B, and C distinguish among the 
different critical densities, and finally the last digit refers to the amplitude of the initial perturbation. 



Model 


p„-/, (gem ^) 


Amplitude a 


Final outcome 


G6A1 


5 X 10-'2 


0.1 


Binary 


G6B1 


5x 10-'* 


0.1 


Quadruple 


G6C1 


5 X 10-'5 


0.1 


Single 


G6B5 


5x lO-" 


0.5 


Single 



The aim of the present models is to explore, using double precision variables, 
the sensibility of fragmentation to both the effect of thermal retardation due to non- 
isothermal heating that is controlled by the critical density, and the effect to increase 
of perturbation amplitude. 

The construction of the initial conditions ensures that the collapse initial stage 
is similar for all models considered here. The initial phase of the collapse proceeds 
in such a way that the material falls down to the rotation plane, while material near 
the central midplane undergoes a weak expansion perpendicular to the rotation axis, 
causing the formation of two overdense blobs from the initial m = 2 perturbation 
seed. At about the end of the first free-fall time, the expansion stops and the middle 
region begins to collapse, the blobs fall toward the center and merge to form a prolate 
structure. By this time, the overall cloud has already been compressed into a flat disk 
with an inner bar that begins to rotate. 
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We begin the discussion with models G6A1, G6B1, seen in figures [T] and |2] The 
m odel G6A1 collapses fas ter than the corresponding single precision model reported 



Arreaga-Garcfal (12007b . but the overall dynamics looks similar. We now proceed 



to decrease the critical density, as shown in model G6B 1 . We find however impor- 
tant differences in comparison to I Arreaga-Garcial (120071 ), which reports as an end 
product a binary system formed by a transient quadrupole system, but we obtain as 
a final product a stable quadrupole system. The fact that by diminishing the criti- 
cal density the f ragmentation enhances , found for Gaussian| Arreaga-Garcfa (2007i); 
Arreaga ( 2008 1 and Plummer models Arreaga-Garcia ( 201oi , seems to be repro- 
duced for models G6A1 and G6B1. However, by diminishing even more the critical 
density in model G6C1, see table [1] the collapse and fragmentation is slowed down 
and the fragmentation is less favored, as seen in figures |2]and [3] 

The model G6B5, shown in figure|4] has a bigger initial perturbation amplitude 
that provokes that systems collapse earlier to form a filamentary structure, which 
evolves to a transient binary system that later collapses in a single system. The 
final system ejects significant amount of gas out of the core. In comparison the 
same model but with a smaller amplitude, model G6B1, forms a stable quadrupole 
system. 



.D6aE-10 gr/cm^ t -1.294at^^ p - 1.525E-09 gr/cm^ t -l.aiarit^^ p - 2.a22E-09 gr/cm" t -l.aB24t^^ 

□ .030 ^^^^^^^^^^^^^^^^^^^^^^H 




i= ISO Log(p'p) 



-o.oao -o.oiQ 0.000 
f= SIC; Log(/; p ) 



j= 475 U^g(p p ) 



Fig. 1 Iso-density contour plots at the equatorial plane of model G6A1 for three different times. 




Fig. 2 Iso-density contour plots at the equatorial plane of model G6B 1 for three different times. 
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13;jE-iy gi/oin-' t -1.7E>01t p - 1.791E-1;;; gr/uiir' t -■.ISnVi.t p - 2.7mE-V,i gr/cm'' t -:;.;il7r3t,. 



-C,04 -0.02 O.D 

N= iOl Log(p/p ) 



-3. 50 0.54 



k= 1D1 log(/)/o ) 



-0.2 -U.i 
k= b?5 J-og{p/p) 



Fig. 3 Iso-density contour plots at the equatorial plane of model G6C1 for three different times. 



D.280E-12 gr/cm^ t -1.5676t 



• 



^.9;^7E-11 gr/cm^ ^1.8244t,, 



Fig. 4 Iso-density contour plots at the equatorial plane of model G6B5 for three different times. 



5 Conclusions 

In this work, we have followed the early phases of cloud collapse and fragmentation 
up to the formation of the proto stellar core using the GADGET-2 code with high 
spatial resolution and double precision, using 10^ SPH particles. The initial condi- 
tions for the cloud model s are cho s en to be the standard isothermal test case, as in 
the variant considered in I Burkerd (|l993l) but for a centrally condensed, Gaussian 
cloud, that was first treated in l Bossl ( ll99lh , and further considered by other authors 
QVrreaga-Garci'a (2007). The main results are summarized as follows: 

By augmenting from single to double precision, the collapse happens earlier and 
the number of e nd products a ugments, as seen in our simulations in comparison the 
same models in l Arreaga-Gar cfa (2007). 

On the other hand, we find that the effect of diminishing the critical density of the 
barotropic equation of state, provokes the collapse to slow down, and this enhances 
the fragments' change to survive. However, this effect happens up to a threshold 
density, as seen in our simulation G6C1, where considered a low critical density, 
Pc = 5 X 10^'^g/cm^, and a single system was formed. 

Moreover, models with a bigger initial perturbation amplitude provoke that sys- 
tems collapse earlier, and in some cases, in the form of a filamentary structure which 
evolves to a transient binary system that later collapses in a single system. 
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